State diagrams for harmonically trapped bosons in optical lattices 
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We use quantum Monte Carlo simulations to obtain zero-temperature state diagrams for strongly 
correlated lattice bosons in one and two dimensions under the influence of a harmonic conflning 
potential. Since harmonic traps generate a coexistence of superfluid and Mott insulating domains, 
we use local quantities such as the quantum fluctuations of the density and a local compressibility 
to identify the phases present in the inhomogeneous density profiles. We emphasize the use of the 
"characteristic density" to produce a state diagram that is relevant to experimental optical lattice 
systems, regardless of the number of bosons or trap curvature and of the validity of the local-density 
approximation. We show that the critical value of U/t at which Mott insulating domains appear 
in the trap depends on the filling in the system, and it is in general greater than the value in the 
homogeneous system. Recent experimental results by Spielman et al. [Phys. Rev. Lett. 100, 
120402 (2008)] are analyzed in the context of our two-dimensional state diagram, and shown to 
exhibit a value for the critical point in good agreement with simulations. We also study the effects 
of finite, but low (T < t/2), temperatures. We find that in two dimensions they have little influence 
on our zero-temperature results, while their effect is more pronounced in one dimension. 
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I. INTRODUCTION 

A great amount of experimental and theoretical 
work [l| has followed the successful realization of the 
superfluid(SF)-to-Mott-insulator transition in ultracold 
bosonic gases trapped in optical lattices, in three 
two Q, and one [j] dimensions. Ultracold atoms on op- 
tical lattices are envisioned as ideal analog simulators 
of Hamiltonians such as the fermion Hubbard model, 
which so far suffers from the lack of reliable analyti- 
cal and numerical results. As a first step toward this 
goal, intensive efforts are currently under way to validate 
this approach by comparing experimental and theoreti- 
cal results for systems such as the Bose-Hubbard model 
that are amenable to both treatments. However, even 
this is challenging, since the phase diagram of the Bose- 
Hubbard model is complicated by issues such as spatial 
inhomogeneities (recognized early in Refs. 0,01), finite- 
temperature effects 0,11, [l^ , and the limited set of ex- 
perimental tools available to characterize those systems. 

The phase diagram of the homogeneous Bose-Hubbard 
model is known to consist of (i) superfluid phases for all 
incommensurate fillings and arbitrary values of the ra- 
tio between the on-site repulsion and the hopping pa- 
rameter (U/t). The system is also superfluid for com- 
mensurate fillings when U/t is smaller than some critical 
value (U/t)c, which depends on the dimensionality of the 
system and on the (integer) filling, (ii) Mott insulat- 



ing phases are present for commensurate fillings when 
U/t > {U/t)c [nl, [H, [H, [H, [3. These two phases 
have been shown to coexist when a confining potential 
is added to the model [1, i, [H, [13, [11. A feature in 
the trap, which is advantageous from the experimental 
point of view, is that in inhomogeneous systems Mott 
insulating domains appear for a broad range of fillings, 
as opposed to the few commensurate fillings required for 
the translationally invariant system. This feature comes 
at a price, sometimes ignored for the sake of simplicity 
by both experimentalists and theoreticians working on 
these systems. In trapped lattice bosons the critical value 
[U /ty^ for the formation of Mott insulating domains in 
different places in the trap not only depends on the local 
filling and the dimensionality of the system (the case for 
homogeneous systems) but also on the total filling N in 
the trap and the curvature of the confining potential V . 
A useful, but approximate, understanding of the effect 
of confinement, which we will not employ here, can be 
obtained through the use of the local-density approxima- 
tion (LDA). 

One may think that having to deal with different fill- 
ings and trapping potential parameters in different ex- 
periments makes the determination of a state diagram 
[l9| in one particular experimental setup not relevant to 
any other. This is not the case. In Refs. [23, [sij, it was 
shown by means of quantum Monte Carlo (QMC) simu- 
lations that for lattice fcrmions in one dimension one can 
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define a scaled dimensionless variable, the characteristic 
density p = NaiV/tY/"^ (where a is the lattice spacing), 
which allows one to build a state diagram in the plane 
p vs U /t that is insensitive to the individual values of 
N and V [l^]- The form for p can most simply be un- 
derstood as arising from dimensional arguments: the trap 
curvature V has units of energy/length^, so that (F/t)^/^ 
has units of inverse length. One can then define a length 
scale ^ = which for trapped systems plays 

a role similar to the system size L in the homogeneous 
case. The characteristic density p = Nhu/^ is then a di- 
mensionless quantity, which for trapped systems is the 
analog of the filling per site n = Ni,a/L in the homo- 
geneous case. In other words, p defines how one should 
approach the thermodynamic limit in trapped systems. 

The characteristic density is not only relevant to one- 
dimensional systems. It can be g eneralized to higher di- 
mensions d, p = Na'^{V/t)'^^^ lUH^I- recent work, the 
state diagram of the three-dimensional Fermi-Hubbard 
model in a harmonic trap was obtained by using a com- 
bination of dynamical mean-field theory (a treatment 
which ignores the momentum dependence of the self- 
energy while retaining its time fluctuations) and the LDA 

One may construct the state diagram in a trap using 
results for the homogeneous system combined with the 
LDA. However, this approach is only approximate for 
finite systems. While LDA has been shown to give rea- 
sonably accurate results in many regimes under a con- 
fining potential, it is certainly bound to fail close to the 
critical values of U /t at which Mott insulating domains 
are formed. This is because in the homogeneous sys- 
tem, there are diverging correlations when one crosses the 
quantum critical region, i.e., finite-size effects in the trap 
not only become relevant but are also unavoidable. We 
will show in this paper that in those regimes where LDA 
fails, our exact QMC-based state diagram is an accurate 
tool for characterization of the experimental results. 

Since state diagrams are not available for lattice 
bosons, in this paper we use world-line quantum Monte 
Carlo simulations to generate the zero-temperature 
state diagram for the one- and two-dimensional (2D) 
Bose-Hubbard model under the influence of a harmonic 
confining potential. The state diagram in two dimen- 
sions allows us to analyze recent experimental results 
for confined two-dimensional bosons in optical lattices 
Q without the need to perform simulations for the 
same parameters as in the experiments. We find that 
the critical values for the formation of Mott insulating 
domains reported in Rcf. 0] are consistent (within 
experimental errors) with our theoretical results for 
inhomogeneous systems. Finally, we will discuss the 
effect that a small increase in the temperature (T < t/2) 
has on our ground-state results. The latter is found to 
have little consequence in two dimensions, where the 
U/t values which dcmark the boundaries between states 
only change by a few percent. 



In one dimension, we find that the trapping potential 
has an even stronger infiuence on the critical values of U jt 
at which Mott insulating domains appear in the trap. We 
will also show that low temperatures have a pronounced 
effect on the inhomogeneous states in the trap. The im- 
portant issue of how to detect experimentally the differ- 
ent phases in the trap is addressed in several recent works 
(see, e.g., [l7, 18, .2|i, i31|, i32|, i33|, iSSij iM]) and will 
not be discussed here. 



II. TRAPPED BOSONS IN TWO DIMENSIONS 

A. Hamiltonian and local observables 

For deep enough optical lattices and low temperatures, 
confined bosons in two dimensions can be described by 
the Bose-Hubbard Hamiltonian [l^ 

H = (ala^.+H.c.) +^^n,(n,-l) 

+Vj2r!n,. (1) 

i 

Here al and are the creation and annihilation oper- 
ators of a boson at site i, located at a distance = 
\Jx\ +y1 from the center of the trap. Xi and yi are 
given in units of the lattice spacing, set to unity in this 
work. Ui = aja^ is the particle number operator. The on- 
site interaction parameter is denoted by C/ (J7 > 0), the 
nearest-neighbor hopping amplitude is denoted by i, and 
V is the curvature of the harmonic confining potential. In 
experiments on optical lattices, [/, and V can be con- 
trolled by changing the intensity of the laser beams that 
produce the lattice. One can also control U separately 
using Feshbach resonances We performed QMC sim- 
ulations in the canonical ensemble, using the world-line 
algorithm [sj . For our-zero temperature study, we have 
taken the inverse temperature to be /3t = 18, which is 
sufficient to obtain ground-state results for the observ- 
ables considered here (see discussion in Fig. [3]) , and for 
discretizing imaginary time we have chosen t At = 0.1. 

In order to create the state diagram, we will moni- 
tor three local observables. The first one is the density. 
Plateaus with constant integer density are the so-called 
Mott insulating domains 5]. However, since the local 
density always crosses integer fillings as the total num- 
ber of bosons in the trap is increased, and since close to 
integer local fillings and for large enough values of U jt 
one always sees some kind of shoulder appearing in the 
density profiles (see, e.g.., Figs. [T] and [3]), identifying the 
formation of Mott domains only by means of the density 
is not accurate [l^, [2lj . 

The second local quantity we monitor is the quantum 
fiuctuations of the density, also referred to as the variance 
of the density, 

A, = {n\) - {n,)\ (2) 
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As shown in Refs. [20|, |2l|, this quantity exhibits a local 
minimum for densities closest to integer fillings. In ad- 
dition, once a Mott insulating domain forms, the value 
at the minimum equals the value of the variance in the 
Mott insulating phase of an equivalent homogeneous sys- 
tem, i.e., a translationally invariant system with exactly 
the same density and U /t [23,[2l|. We will show here that 
in our finite lattice-boson systems (such as the ones cre- 
ated experimentally), shoulders with ~ I can appear 
without the variance of the density on those shoulders 
attaining the value in the homogeneous Mott insulator. 
That is, those shoulders are not local Mott insulating 
domains. 

The final quantity that we measure here is the local 
compressibility defined in Refs. [sl. [l8|. 



(a) 





(3) 

This local compressibility quantifies the response in the 
on-site density to a local change of the chemical poten- 
tial. Ki also exhibits a minimum around integer fillings, 
and it behaves, qualitatively, similarly to A^. An anal- 
ogous distinction between equal-time fluctuations and 
susceptibilities which include unequal-time correlations 
plays an important role in the fermion Hubbard model, 
where the local susceptibility carries additional informa- 
tion about the (Kondo) screening beyond that contained 
in the equal-time local moment. 



B. Density, variance, and local compressibility 

In Fig. [1] we show two density profiles in a two- 
dimensional harmonic trap for values of the on-site re- 
pulsive interaction right before ([//< = 17.5) and right 
after {U/t = 18.5) the formation of the Mott insulator in 
the middle of the trap. For U/t = 17.5 [Fig. [TJa)], one 
can see that a shoulder with density n ~ 1 has devel- 
oped in the system, but the density in the center of the 
trap is shghtly greater than 1. For U/t = 18.5 [Fig.IHb)] 
only a plateau with density n = 1 is seen at the central 
region of the system. Positions in the trap are normal- 
ized by a length scale ^ = \ftfV , which is introduced 
in Hamiltonian ((T|) by the harmonic confining potential. 
With this normalization, systems with different trapping 
potentials and fillings, but with the same characteristic 
density, have identical density profiles [10, [2l|, [2^ . 

Intensity plots of the variance and local compressibil- 
ity profiles corresponding to the systems in Fig. [T] are 
presented in Fig. [2] This figure shows that both A and k 
exhibit minima whenever the density is closest or equal 
to 1. In addition, in Figs.[2]Ja) andl^Jb) one can see that 
the region with rt > 1 at the center of the trap in Fig.[TJa) 
is signaled by local maxima of A and k. 

A more quantitative understanding of the behaviors of 
n. A, and k can be gained by plotting the changes in 
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FIG. 1; (Color online) Two-dimensional density profiles for 
systems with p = 30 (Nt = 1200 and V/t = 0.025), (3 = 18, 
Lx ~ Ly — 60 (number of lattice sites along x and y), and 
in-site repulsion (a) U/t = 17.5 and (b) U/t — 18.5. Positions 
in the trap are normalized by the length scale ^ = \/t/V (see 
text). For the homogeneous model (V/t=0) the critical value 
{U/t)c for the SF-Mott transition in d = 2 is 16.74 [ist . 



these quantities along one spatial dimension, while the 
coordinate in the other spatial dimension is fixed to the 
center of the trap. The results for the systems in Figs. [T] 
and [5] are shown in Fig. [31 which confirms that, indeed, 
minima of A and k correspond to regions with n ~ 1 
in the density profiles. Along with the results in the 
trap, we have plotted as horizontal lines the values of 
the density, variance, and compressibility in the Mott 
insulating phase of the homogeneous system for exactly 
the same values of U/t. 

Figures[3l^a)--(c) {U/t = 17.5) show an important prop- 
erty of these finite trapped systems. While the density 
profile can exhibit a shoulder with n ~ 1 [Fig.[31[a)] when 
crossing n = 1, the variance and local compressibilities 
on this shoulder do not reach their values in the first lobe 
of homogeneous systems, for exactly the same values of 
U /t. The LDA is not valid in this region of the trap. 
Notice that in Figs. [2Ja)-[3]Jc) such a shoulder is present 
for a value of U/t that is greater than the critical value 
{U /t)c for the formation of the Mott insulator in the ho- 
mogeneous system, which is {U /t)c = 16.74 for n = 1 in 
two dimensions [15j. Therefore, no local Mott insulating 
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FIG. 2: (Color online) Intensity plots of the local quantum 
fluctuations of the [(a) and (c)] density A and [(b) and(d)] 
local compressibility k for [(a) and (b)] U/t = 17.5 and [(c) 
and (d)] U /t = 18.5. In these systems p = 30 (iVb = 1200 and 
V/t = 0.025), fi = 18, and = Ly ^ 60. The corresponding 
density profiles are plotted in Fig. [T] 



phase is found in a trap for U/t > [U /t)c and the density 
at the center of the trap is n > 1. In an experiment with 
bosons trapped in optical lattices, the critical value for 
the formation of a Mott insulating state may be greater 
than the critical value in homogeneous systems. This is 
due to the finite curvature of the trapping potential and 
the finite (and sometimes small) extent of the system. 
How much the critical value is shifted from the homoge- 
neous prediction is something that, as we will show later, 
will strongly depend on the dimensionality of the system. 

In Figs.[i;d)-[3i;f) one can see that for U/t ^ 18.5 a full 
region with n = 1 occupies the center of the trap. In that 
region, A and k have exactly the same values as that in 
the homogeneous system with n = 1 for an identical on- 
site repulsion U jt. Hence, we conclude that the domain 
in the middle of the trap is in this case a Mott insulator. 

To conclude this subsection on local observables, we 
discuss whether the lattice size considered in Figs. [T] and 
[5]is big enough so that boundary conditions are irrelevant, 
and whether the temperature is low enough so that n. A, 
and K behave as they will in the ground state. In Fig. [31 
we present results for a system with identical Hamilto- 
nian parameters and filling but at a lower temperature 
/3 = 24 and in a smaller lattice with = Ly = 50 sites in 
each direction. They can be seen to be essentially indis- 
tinguishable from those with /? = 18 and — Ly — 60. 
Hence, at least for our local observables, the tempera- 
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FIG. 3: (Color online) Cuts across the center of the trap for 
the two-dimensional density, variance, and local compressibil- 
ity plots in Figs.[T]and[2][p = 30 (Nt = 1200 and V/t = 0.025), 

l3 = 18, and L X Ly = 60]. Horizontal lines show the results 

in homogeneous systems for the same values of U/t, jS, and 
n = 1. These results emphasize that the LDA does not hold 
in some regions of the trap. In the n ^ 1 domains in panel 
(a) that develop around commensurate filling, the variance in 
panel (b) does not take on the value predicted by the homo- 
geneous system for n = 1. We also plot results for a smaller 
system with = Ly = 50, at lower temperature /3 — 24, and 
the same value of p = 30 [Nb = 1200 and V/t = 0.025). 



ture considered is low enough and the systems sizes are 
big enough so that they have no influence in our results. 



C. State diagram 

Uniform systems of different sizes have identical prop- 
erties as long as they are chosen to have the same den- 
sity p — Nb/L'^, or in a lattice the same n = N^a'^/L'^, 
and are sufficiently large. In a similar way, the charac- 
teristic density generalizes this for confined systems and 
allows one to obtain state diagrams that are independent 
of particular choices of boson number and trap curvature 
provided the systems are also large. Finite systems, such 
as the ones realized experimentally, pose an additional 
complication as finite-size effects become relevant. They 
are most important close to where local phases (such as 
those described in Sec. Ill Bp appear for the first time in 
the trap. This implies that, as also shown in Sec. Ill Bl 
the LDA may fail to describe the actual density profiles 
in an ultracold gas experiment on an optical lattice. 

In this section we show that, for systems like the ones 
realized experimentally, a state diagram in the plane p 
vs U /t can accurately predict the local phases in trapped 
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experiments even when the LDA fails. What this means 
is that for those system sizes (where the occupied part 
of the lattice has 50 < L^, Ly < 100), changing the spe- 
cific trap parameters does not appreciably change the 
observed local phases under the confining potential. In 
order to create such a state diagram, we have considered 
systems with two different values of the trap curvature 
and many different fillings. We fixed the inverse tem- 
perature to /3 = 18 and the largest lattices considered 
here had = Ly ~ 64 lattice sites. The lattice sizes 
considered here are comparable to the ones realized ex- 
perimentally 3. 

In the left panels of Fig. Il[(a)-(c)], we compare two 
systems with the same characteristic density p = 20, but 
two different trapping potentials and fillings [iVf, = 800, 
V/t = 0.025 and Nb = 500, V/t = 0.04]. Our results 
for the scaled density, variance, and local compressibil- 
ity profiles are essentially indistinguishable between those 
two traps, despite the large difference in particle number 
and curvature. In Figs. [Ha)-l4jc), we have chosen one 
of the lowest characteristic densities that support a Mott 
insulating state in the trap. With increasing p, at the 
same values of V/t, the differences between the scaled 
profiles in the two traps become even smaller. 

In the right panels of Fig. |4[(d)-(f)], we compare two 
systems in which an additional superfluid phase with 
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FIG. 4: (Color online) Density, variance, and local compress- 
ibility in one direction across the center of the trap. We com- 
pare systems with the same characteristic density but dif- 
ferent trapping potentials and fillings. Left panels [(a)-(c)]: 
U/t = 24, p = 20, with Nb = 800, V/t = 0.025 and Nb = 500, 
V/t = 0.04. Right panels [(d)-(f)]: U/t = 19, p = 35, with 
Nb = 1400, V/t = 0.025 and Nb = 875, V/t = 0.04. In all 
cases P = 18. Horizontal lines show the results for homoge- 
neous systems with n = 1 and the same values of U/t and (3 
as in the trap. 



n > 1 is present inside a Mott insulating domain with 
n = 1. This state has a richer structure than the one in 
the left panel. The characteristic density is the same in 
both systems, p = 35, but the trap curvature and fillings 
are different [Nb = 1400, V/t = 0.025 and TV^ = 875, 
V/t = 0.04). We have chosen the ratio of U /t in this case 
to be very close to the critical value for the formation 
of the Mott insulator at the center of the trap. Figures 
llld)-[Hf) show that even in this more complicated case, 
with several different domains and close to a transition 
between states, the scaled results for our local observables 
are also almost indiscernible even for two very different 
systems, as long as they have the same characteristic den- 
sity. 

In Figs. [3]Ja)-131[c), we have shown results for local 
quantities in a system that cannot be described within 
the LDA. In Fig. O we provide an additional example 
where two different confined systems with U/t> {U/t)c 
[{U/t)c from the homogeneous case] fail to exhibit Mott 
insulating domains. This figure also shows that even 
though the LDA is clearly not applicable to those sys- 
tems, the density profiles and local compressibilities in 
both of them, which have the same characteristic den- 
sity, are almost indistinguishable. As in previous figures, 
the number of particles in the largest system is almost 
twice that of the smallest. These results confirm that 
the characteristic density, being the proper quantity to 
define the thermodynamic limit in a trapped system, pro- 
vides genuine guidance for finite systems when the LDA 
is not valid. 

Following the discussion in Sec. Ill Bp . we locate Mott 
insulating domains in the trap by comparing the variance 
of the density and local compressibilities in the confined 
system with those of a Mott insulator in the homogeneous 
case. Whenever they coincide, wc conclude that a local 
Mott insulating domain is present. With this criterion, in 
most cases Mott insulating domains emerge in the trap 
for the same value of U /t independent of whether one 
considers A or k. In a few cases, especially at low fillings, 
we found A in the trap to agree with the value in the 
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FIG. 5: (Color online) (a) Density and (b) local compressibil- 
ity in one direction across the center of the trap. We compare 
systems with the same characteristic density p = 35 but differ- 
ent trapping potentials and fillings of V/t = 0.025, A^;, = 1400 
and V/t — 0.04, Nb = 875. Horizontal lines show the results 
for homogeneous systems with n = 1 and the same values of 
U/t and /3 = 18 as in the trap. 
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FIG. 6: (Color online) State diagram for bosons in a two- 
dimensional lattice in a harmonic confining potential. The 
boundaries between states were determined using two differ- 
ent traps, with V/t = 0.025 (filled symbols) and V/t = 0.04 
(open symbols). Both traps lead to very similar results except 
for some small differences at the lowest characteristic densi- 
ties. The states are; (I) a pure superfluid phase [Figs. EJa)- 
|3jc)], (11) a Mott insulating phase at the center of the trap 
surrounded by a superfluid phase with n < 1 [Figs. |3Id)-[3If ) 
and Figs. |4][a)-[4{c)], (III) a superfluid phase with n > 1 at 
the center of the trap surrounded by a Mott insulating phase 
with n = 1, and an outermost superfluid phase with n < 1 
[Figs. |4jd)-|4jf)]. The vertical dashed line signals the criti- 
cal value of U/t for the formation of the Mott insulator with 
n = 1 in the homogeneous case [l^ . The inverse temperature 
in the trapped systems is /3 = 18. 

homogeneous case before k did, usually by a difference 
S{U/t) ~ 0.5. In those cases, the value of U/t reported 
at the state boundary is the average between the one 
obtained by A and the one obtained by k. Note that this 
discrepancy represents, at worst, an uncertainty in the 
boundary of only a few percent. 

Figurc[n]is the central result of this paper: the state di- 
agram for lattice bosons in a two-dimensional lattice un- 
der the influence of a harmonic trap. The different states 
depicted represent the following. (I) Only a superfluid 
phase is present in the trap, the situation in Figs. [2{a)- 
^c). (Notice that in this state the density at the center 
of the trap can take arbitrarily high values.) (II) A Mott 
insulating phase at the center of the trap is surrounded 
by a superfluid phase with n < 1. Examples of this state 
are given in Figs. [3i;d)-[3i;f) and Figs. [4i;a)-l4i;c). (Ill) A 
local superfluid phase with n > 1 is present at the center 
of the trap. This phase is surrounded by a Mott insu- 
lating domain with n = 1, inside a superfluid phase with 
n < 1. An example of that state can be found in the 
right panel of Fig. |li:d)-[l);f). 

An important feature of the state diagram of Fig. [HI 
is that except for a small window of characteristic den- 
sities (p ~ 23), where the critical value of U/t for the 
formation of the Mott insulating domain at the center 
of the trap {U/t)]. is very close to the value in homo- 



geneous systems {U/t)c, for most characteristic densities 
{U/t)l departs significantly from {U/t)c] i.e., this state 
boundary depends strongly on p. At least for fermions 
[sst . this is the state boundary that currently can be 
accurately determined experimentally. This is because 
the total double occupancy, which is a good indicator for 
distinguishing states with a Mott insulator at the center 
of the trap from states with higher density [1^, can be 
accurately measured in such systems [ssj . Similar mea- 
surements in bosonic systems, which have site occupancy 
higher than two, could allow experimentalists to repro- 
duce our results over an even broader range of the state 
diagram. 

The other state boundary of interest is the one be- 
tween states I and III. This boundary, on the other hand, 
does not depend strongly on the characteristic density 
[3^. This means that by experimentally measuring the 
boundary between states I and III, which is sensitive to 
transport and coherence probes, one would obtain a crit- 
ical value of U/t for the formation of the Mott insulator 
that is close to the one in the homogeneous system. 

To conclude, we should add that with increasing char- 
acteristic density new states with Mott insulating do- 
mains with n > 1 appear in the trap. Those states will 
form for larger values of the in-site repulsive interaction, 
and will be more difficult to detect experimentally. 



D. Comparison with NIST experimental results 

As mentioned in Sec. HI recent experiments have 
achieved the superfiuid-to-Mott-insulator transition in 
two dimensions [3| . Here we provide a quantitative com- 
parison with our state diagram. 

We use the precise values of the curvature V gener- 
ated by the optical and magnetic traps, and interaction 
and tunneling energies U and t (computed using a band 
structure calculation) appropriate to the NIST experi- 
ments . The final parameter of interest is the filling. 
The NIST experiments employed an array of about 70 
independent 2D systems of different filling. The number 
of particles reported for the central 2D slice (the one with 
maximal filling) Q was Ni, w 3000, while the average was 
Nb ~ 2000. 

Figure [7] shows the trajectories of systems with differ- 
ent fillings, and the NIST trapping and energy parame- 
ters, in the QMC state diagram. We have also plotted 
the experimentally reported critical values of U/t for the 
formation of the Mott insulator obtained using the con- 
densate fraction and the FWHM of the momentum dis- 
tribution function. For intermediate fillings (Ni, between 
1000 and 2000 particles) corresponding to the average 
filling of the array of 2D layers, the experimental results 
intercept (within their error bars) the QMC boundary 
between states I and II. The error bars reported in the 
experiments arise mainly from the uncertainty in the lat- 
tice depth, which translates into around a ±10% uncer- 
tainty in the ratio U/t. For the largest fillings of the 2D 
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FIG. 7: (Color online) State diagram in Fig. [6] with exper- 
imental trajectories for systems with different fillings, A*'(, = 
1000, 2000, and 3000. V, t, and U have been chosen to cor- 
respond to the specific experiments is Ref. 0|. Vertical lines 
depict the critical values oi U/t reported in Ref. computed 
using the condensate fraction in one case and the full width at 
half maximum (FWHM) of the momentum distribution func- 
tion in the other [3]. The experimental errors are also plotted. 
The experiment and QMC give values for the critical coupling 
in excellent agreement. 



cter, so we consider here temperatures T < t/2. The 
effect of higher temperatures {t < T < U) on. the Mott 
insulating state in trapped systems has been discussed in 
previous works (see, e.g., Refs. 0,[i,i,[i3). 

Our QMC simulations for two-dimensional systems 
with T < t/2 show that low temperatures have a small 
effect on the state diagram in Fig. [51 We have found that 
in general the state boundaries move to larger values of 
U /t, with a change of up to 4% of the ground-state value 
depending on the filling. This effect is particularly ev- 
ident for the lowest characteristic densities. For larger 
characteristic densities, the ones that support state III 
in Fig. [HI we find the effect of low temperatures to be 
less important. 

As an example of our results, in Fig. [5] we compare 
density profiles and compressibilities in confined systems 
with p = 44 at different temperatures and interaction 
strengths. One can see there that while density and com- 
pressibility profiles in state I are almost indistinguishable 
for the temperatures considered here [Figs.[5]Ja) and[5]^b), 
where U/t = 16], the ones in state III can exhibit ap- 
parent changes when close to the boundary with state II 
[Figs.[5Jc) and[Hfd), where U/t = 21]. However, with just 
a small increment of the interaction strength so that one 
creates state II in the ground state [Figs. [H^e) and[8{f). 



slices {Ni, between 2000 and 3000 particles), we find that 
the experimentally reported critical values are very close 
to the boundary between states I and III. 

Overall, the agreement between our numerical calcula- 
tions and the experimental results is remarkable. It could 
even be considered surprising, taking into account that 
in the experiments only global quantities such as the con- 
densate fraction and the FWHM of the momentum dis- 
tribution function have been measured, while our state 
diagram is based on the behavior of local quantities such 
as the variance of the density and the local compress- 
ibility. The agreement between those two approaches is 
a priori not guaranteed [l^. New experiments in which 
the filling in the 2D sections is better controlled could al- 
low experimentalists to map the state diagram in Fig. [6] 
Given the presence of the trapping potential in the ex- 
periments and the fact LDA fails close to {U/t)c of the 
homogeneous system, we believe this is a better strat- 
egy to follow to validate experiments than trying to map 
the phase diagram of the homogeneous system. Recent 
progress on the direct imaging of the density profile will 
also allow further detailed comparisons with simulation. 



E. Finite temperature 

In order to better compare with the experiments, it 
is important to understand how a small increase in the 
temperature affects the zero-temperature state diagram 
reported in this paper. By small we mean a temperature 
that is low compared to the value of the hopping param- 
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FIG. 8: (Color online) Comparison between local quantities 
in systems with different temperatures. We have plotted the 
density (left panels) and local compressibility (right panels) in 
one direction across the center of the trap. The ratio between 
the on-site repulsion and the hopping parameter is increas- 
ing from top to bottom: [(a) and {h)]U/t = 16, [(c) and (d)] 
U/t = 21, and [(e) and (f)] U/t = 21.5. All systems have the 
same characteristic density p = 44 [Ni, = 1100, V/t = 0.04). 
Horizontal lines show the results for homogeneous systems 
with n — 1, and the same values of U/t as in the trap. The 
differences bewteen k, in homogeneous systems with I3t = 2, 6, 
and 18 are indistinguishable in the figure. 
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where U/t = 21.5] the system has a full Mott insulating 
domain at the center of the trap for /3 = 18 and /3 = 6, 
while (3 — 2 \s very close to forming one. The generic ef- 
fect of the temperature can be seen to be a displacement 
of the state boundaries to larger values oi U/t. 



III. TRAPPED BOSONS IN ONE DIMENSION 

In one dimension, previous work has already shown 
that systems with the same characteristic density p = 
Nb \/V/t have identical re-scaled density profiles [23 , so 
here we will not present a detailed discussion (i.e., the 
d = 1 analogs of Figs. [31 [H and [5] presented for d — 2). In 
addition, the behavior of the local observables defined in 
Sect.[lllis qualitatively similar to that in two-dimensional 
systems. Therefore, in this section, wc use the same cri- 
teria to create the d = 1 state diagram, depicted in Fig.[9l 
The labeling of the states follows the same convention as 
in two dimensions. 
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FIG. 9; (Color online) State diagram for bosons in a one- 
dimensional lattice with a harmonic confining potential. The 
state diagram was built using a trap with V jt = 0.008, L = 
100, and /3 = 10. The states are designated in the same way as 
in the two-dimensional case: (I) a pure superfluid phase, (II) 
a Mott insulating phase at the center of the trap surrounded 
by a superfluid phase with n < 1, (III) a superfluid phase 
with n > 1 at the center of the trap surrounded by a Mott 
insulating phase with n = 1, and an outermost superfluid 
phase with n < 1. The vertical dashed line signals the critical 
value of U/t for the formation of the Mott insulator with n = 1 
in the homogeneous case [T^ . 

Figure [21 shows that in one dimension the trapping po- 
tential has a more pronounced effect on displacing the 
critical value for the formation of the Mott insulator to- 
ward larger values of U/t. This can be understood if one 
compares the homogeneous Mott lobes in one dimension 
with those of higher dimensional systems, and is some- 
how similar to what was found in Refs. [lO, Hlj for the 
fermionic case, where, in the trap, the lowest value of 
U/t at which a local Mott insulating phase was found 
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FIG. 10: (Color online) Comparison between local quantities 
[density (left panels) and local compressibility (right panels)] 
in one-dimensional systems with different temperatures. The 
ratio between the on-site repulsion and the hopping parameter 
is increasing from top to bottom: [(a) and (b)] U/t = 10, [(c) 
and (d)] U/t = 14, and [(e), (f)] U/t = 16. AU systems have 
the same characteristic density p — 6.96 {Nt = 55, V/t = 
0.016). Horizontal lines show the results for homogeneous 
systems with n = 1, and the same values of U/t as in the 
trap. The differences between k in homogeneous systems with 
pt = 2, 4, 6, and 10 are indistinguishable in the figure. 



was U/t 3, to be compared with {U/t)c = in the ho- 
mogeneous case. Another important difference between 
the state diagrams in two and one dimensions is that in 
the latter all states boundaries have a stronger depen- 
dence on the characteristic density. 

In order to study the effect of finite temperatures in 
one-dinicnsional systems, we performed quantum Monte 
Carlo simulations using the stochastic Green's-function 
algorithm [4l[. This algorithm allows us to compute the 
momentum distribution function of the bosons, a quan- 
tity that is very difficult to compute with the world-line 
approach. 

In Fig. [TOl we show examples of the effect of the tem- 
perature in one dimension. The temperature is increased 
up to T = t/2, similarly as in the two-dimensional case. 
Figures flW a) and fTHT b) depict a system that at zero tem- 
perature is in state III in Fig. [9l Increasing the tempera- 
ture from /3 = 10 to /3 = 2 in that state does not produce 
large changes in the density and compressibility profiles. 
The region with n ~ 1 is the most affected by the tem- 
perature, as the Mott plateaus to the side of the central 
region with n > 1 have shrunken. 

Figures [TOFc) and [TOFf) deal with the state that at 
zero temperature has a Mott insulator at the center of 
the trap (state II in Fig. For U/t = 14 [Figs. IMc) 
and lior d)]. one can see that increasing the temperature 
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beyond /3 = 6 melts the Mott insulator at the center of 
the trap, producing a central region with n > 1. As one 
increases U/t, U/t ^ 16 [Figs. fTOT c) and flOT f)]. one can 
see that the central Mott insulating domain survives for 
all temperatures analyzed here. 

Overall, Fig. [10] shows that the effect of low temper- 
atures in the state diagram in one dimension is qualita- 
tively similar to that in two dimensions; i.e., the bound- 
aries of state II move to larger values of U/t. How- 
ever, our simulations show that one-dimensional systems 
are more affected by finite temperatures than their two- 
dimensional counterparts, so that any experimental com- 
parison with our state diagram in Fig. [S] would require 
very low temperatures. These results also mean that the 
critical values for the formation of the insulator in one- 
dimensional trapped systems at finite temperatures arc 
further away from the critical value in the homogeneous 
system at zero temperature. 

To conclude, we show in Fig. [Tl] the effect of the tem- 
perature on the momentum distribution function of the 
systems with U/t = 10 and U/t = 14 in Fig. [lOl (The 
momentum distribution function for the systems with 
U/t = 16 remains almost unchanged for the tempera- 
tures considered here.) Figure [TlT a) shows that when 
most of the system is in a superfluid phase at zero tem- 
perature, increasing the temperature leads to a reduction 
in the zero-momentum peak occupation, which is a con- 
sequence of the reduction in the correlation length in the 
superfluid domains. On the other hand, Fig.fTlTb) shows 
that when the ground state of the system has a very large 
Mott insulating domain at the center of the trap, an in- 
crease in the temperature can increase the occupation of 
the zero-momentum peak when the Mott insulator melts 
to give rise to a region with n > 1 at the center of the 
trap. 
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FIG. 11: (Color online) Momentum distribution function in 
one- dimensional systems with different temperatures. The 
ratio between the on-site repulsion and the hopping parameter 
are (a) U/t = 1G (at zero temperature the system is in state 
III), and (b) U/t = 14 (at zero temperature the system is 
in state II). All systems have the same characteristic density 
p = 6.96 {Nb = 55, V/t = 0.016). 



IV. CONCLUSIONS 

In this paper we have computed the state diagrams 
for the formation of Mott insulating domains with n = 
1 in lattice-boson systems confined by harmonic traps. 
Results have been presented in one and two dimensions 
and the lattice sizes considered had between 50 and 100 
sites in each direction, values which are comparable to 
those studied by the NIST group. 

We have shown that for the system sizes considered, 
state boundaries are accurately determined by the char- 
acteristic density, providing a state diagram for different 
experimental choices of trap curvature and particle num- 
ber. Key features of these state diagrams are that the 
critical values of U /t for the formation of a Mott insu- 
lator in a trap measured experimentally need not be the 
same as in homogeneous systems and that the LDA does 
not hold close to the state boundaries. We find that in 
two dimensions, the lowest value of U/t that supports a 
local insulator in the trap is U/t = 17.4, which is ap- 
proximately 4% greater than the critical value in the ho- 
mogeneous case, {U/t)c = 16.74. We have also shown 
that {U/t)'^ for the formation of the Mott insulator at 
the center of the trap increases with the characteristic 
density. These results could be verified experimentally 
by measuring double or higher occupancy in bosonic sys- 
tems. 

A comparison of our two-dimensional results with the 
ones of experiments at NIST shows that the particular 
trap, lattice depth, and particle filling used there give 
trajectories in the state diagram which intersect the crit- 
ical coupling line near its extremal value of U/t. The 
transition point reported is in good agreement with our 
QMC simulations, which explicitly include the confining 
potential. We have also shown that, in two dimensions, 
finite (but low, T < t/2) temperatures have little effect 
on our ground-state results. 

In one dimension, we find that the trapping potential 
has a stronger effect in moving the critical values for the 
formation of the Mott insulator toward larger values of 
U/t. Actually, the lowest value of U/t at which we find 
a local Mott insulating phase in the trap (with up to 
100 sites) is U/t = 5.5, that is, a 52% increase over the 
critical value in the homogeneous case, {U/t)c = 3.61. 
We also find that state boundaries in one dimension are 
more sensitive to the characteristic density in the trap 
and to low temperatures than those in two-dimensional 
systems. 
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